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Abstract 


A third-order .Energy Stable Weighted Essentially Aon- Oscillatory (ESWENO) finite difference scheme 
developed by Yamaleev and Carpenter (AIAA 2008-2876, 2008) was proven to be stable in the energy 
norm for both continuous and discontinuous solutions of systems of linear hyperbolic equations. Herein, a 
systematic approach is presented that enables “energy stable” modifications for existing WENO schemes 
of any order. The technique is demonstrated by developing a one-parameter family of fifth-order upwind- 
biased ESWENO schemes; ESWENO schemes up to eighth order are presented in the appendix. New 
weight functions are also developed that provide (1) formal consistency, (2) much faster convergence for 
smooth solutions with an arbitrary number of vanishing derivatives, and (3) improved resolution near 
strong discontinuities. 
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1 Introduction 

A new, third-order weighted essentially nonoscillatory scheme (called Energy- Stable WENO) has recently 
been proposed and developed by the authors of the present paper [1], In reference [1], we prove that 
the third-order ESWENO scheme is energy stable, that is, stable in an L 2 energy norm, for systems 
of linear hyperbolic equations with both continuous and discontinuous solutions. Stability is explicitly 
achieved (by construction) by requiring that the ESWENO scheme satisfies a nonlinear summation-by- 
parts (SBP) condition at each instant in time. Thus, L 2 strict stability is attained without the need 
for a total variation bounded (TVB) flux reconstruction or a large-time-step constraint [2], [3] and [4]. 
Herein, we generalize and extend the third-order ESWENO methodology [1] to an arbitrary order of 
accuracy. Similar to the third-order ESWENO scheme, the new families of higher order (up to eighth 
order) ESWENO schemes are provably stable in the energy norm and retain the underlying WENO 
characteristics of the background schemes. Numerical experiments demonstrate that the new family 
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of ESWENO schemes provides the design order of accuracy for smooth problems and delivers stable 
essentially nonoscillatory solutions for problems with strong discontinuities. 

Another issue that is addressed in this paper is the consistency of the new class of ESWENO schemes. 
The consistency of any WENO-type scheme fully depends on a proper choice of the weight functions. On 
one hand, for smooth solutions the weights should provide a rapid convergence of the WENO scheme to 
the corresponding underlying linear scheme. On the other hand, the weights should effectively bias the 
stencil away from strong discontinuities. The high-order upwind-biased WENO schemes with conven- 
tional smoothness indicators that are presented in reference [2] are too dissipative for solving problems 
with a large amount of structure in the smooth part of the solution, such as direct numerical simulations 
of turbulence, or aeroacoustics[5], [6]. Furthermore, as has been shown in references [7] and [8], the clas- 
sical weight functions of the fifth-order WENO scheme fail to provide the design order of convergence 
near smooth extrema, where the first derivative of the solution becomes equal to zero. New approaches 
are proposed in references [7] and [8] to improve the error convergence near the critical points. Although 
these new weight functions recover the fifth order of convergence of the WENO scheme near smooth 
extrema, the problem persists if the first- and second-order derivatives vanish simultaneously [8]. An 
attempt to resolve this loss of accuracy is presented in reference [8]. This proposed resolution provides 
only a partial remedy for the problem; the same degeneration in the order of convergence occurs if at 
least the first three derivatives become equal to zero. To fully resolve this problem, we propose new 
weights to provide faster error convergence than those presented in reference [8], and impose some con- 
straints on the weight parameters to guarantee that the WENO and ESWENO schemes are design-order 
accurate for sufficiently smooth solutions with an arbitrary number of vanishing derivatives. 

This paper is organized as follows. In section 2, energy estimates for the continuous and corresponding 
discrete wave equations are presented. In section 3, we present a one-parameter family of fifth-order 
WENO schemes; one value of the parameter yields a central scheme that converges with sixth-order 
accuracy. In section 4, we present a systematic methodology for constructing ESWENO schemes of any 
order and demonstrate the methodology by transforming the family of WENO schemes presented in 
section 3, into a family of fifth-order ESWENO schemes. In section 5, we analyze the consistency of the 
new class of ESWENO schemes and we derive sufficient conditions for the weights functions that ensure 
that the ESWENO schemes are design-order accurate regardless of the number of vanishing derivatives in 
the solution. The tuning parameters in the weight functions are also optimized. In section 6, we present 
numerical experiments that corroborate our theoretical results. We summarize and draw conclusions in 
section 7. 


2 Energy Estimates 

Consider a linear, scalar wave equation 

I y + fy = °- / = au, t > °, ° < £ < i, 

u(0, x) = uq{x) ' 

where a is a constant and uo(x) is a bounded piecewise continuous function. Without loss of generality, 
assume that a > 0, and further assume that the problem is periodic on the interval 0 < x < 1. Applying 
the energy method to equation (1) leads to 


>.«!, = » (2) 

where || • \\l 2 is the continuous L 2 norm. Thus, the continuous problem defined in equation (1) is neutrally 
stable. 

We now develop using mimetic techniques (see [1] or [9]) a class of discrete spatial operators that is 
neutrally stable or dissipative. The continuous target operator used for this development is the following 
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singular perturbed wave equation: 


N 

f + if=E i (- 1 f = au, t> 0, 0 < x < 1, (3) 

m(0,£) = u 0 (x) , 

where fi = fi(u) is a non-negative C°° function of u. As before, we assume that equation (3) is subject to 
periodic boundary conditions. Our goal is to match each spatial term in equation (3) with an equivalent 
discrete term that maintains neutral stability (or dissipates) of the discrete energy norm. 

We begin by showing that the terms on the right-hand side of equation (3) are dissipative thereby 
ensuring stability. Multiplying equation (3) by u and integrating it over the entire domain yields 


1 d_ 

2 dt 


IMIL 


2 aU 


1 N 

= £ 

0 n— 1 ' 


z ,\n — 1 w w u ' 7 

(— 1) u — — fi — — dx . 

' ' 0 ^. 77 , ‘ £»~, 77 . 


dx 


dx r 


( 4 ) 


Integrating each term on the right-hand side by parts and accounting for periodic boundary conditions 
yields the following energy estimate: 


-I 

dP 


N 


= - 2 £ / m(w) 


n— 1 \ 


d n u 

dx n 


dx < 0 


( 5 ) 


All the perturbation terms included in equation (3) provide dissipation of energy. 

Turning now to the discrete case, we define a uniform grid Xj = jAx, j = 0, J, with Ax = 1/J. 
On this grid, we define a flux f = au and its derivative fj, = au s , where u = [u(xo,t), . . . ,u(xj,t)] T 
and u z = [u x (xo,t), . . . , u x (xj,t)] T are projections of the continuous solution and its derivative onto 
the computational grid. Next, we define a pth-order approximation for the first-order derivative term in 
equation (1) as 

Of 

— = Di + 0(Ax p ) . (6) 

ox 

Placing a mild restriction on the generality of the derivative operator (see [1] or [9]), the matrix D can 
be expressed in the following form: 


D = P~ 1 [Q + R } ; Q + Q T = 0 

R=R t ; v t Av > 0 (7) 

P = P T ; v T Pv > 0 

for any real vector v/0. By choosing the matrix R as a discrete analog of the dissipation operator in 
equation (3), we have 

N 

r = J2 d i n M D i n } T , ( 8 ) 

71=1 

where A is a diagonal positive semidefinite matrix and D\ is the difference matrix 




0 

-1 1 

0 


By using the SBP operators [eqs. (6)— (8)], the semi-discrete counterpart of equation (1) becomes 


<9u 

~dt 


+ P~ l Q f 


N 

^p- 1 D 1 "A[Di n ] T f, 

71=1 


( 9 ) 
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where f = au, u = [uo(t), iti(i), . . . , uj(t), } T is the discrete approximation of the solution u of equation 
(1), and Q and A are nonlinear matrices (i.e. , Q = Q( u) and A = A(u). To show that the above 
finite-difference scheme is stable, the energy method is used. Multiplying equation (9) with P yields 


Id 
2 dt 


N 


u||p + au T Qu = -a^2 ([-Di"] T u) A[Di”] t u, 


( 10 ) 


n— 1 


where || • \\p is the P norm (i.e., ||u||p = tr P u). Adding equation (10) to its transpose yields 

7 N 

— ||u||p + au T ( Q + Q t ) u = -2 ([A>i”] t u) T A[P>i n ] T u . 


(ID 


If we account for periodic boundary conditions and the skew-symmetry of Q , then the second term on 
the left-hand side vanishes, and the energy estimate becomes 


N 


dt. 


|u||?, = -2 a]T (m T u) A[Z?rfu < 0 . 


( 12 ) 


71 = 1 


The right-hand side of equation (12) is nonpositive because the diagonal matrix A is positive semidefinite 
[v T Av > 0 for all real v of length ( J + 1)] and a > 0; thus the stability of the finite-difference scheme 
given by equation (9) is assured. This result can be summarized in the following theorem: 

Theorem 1. The approximation [eq. (9)] of the problem [eq. (1)] is stable if equations [(6)-(8)f hold. 

Remark 1. Despite the fact that the initial boundary value problem [eq. (1)] is linear, the finite-difference 
scheme [eq. (9)] constructed for approximation of equation (1) is nonlinear, because the matrices Q and 
A (and in principle P) are assumed to depend on the discrete solution u. 

Remark 2. The only constraints that are imposed on the matrix Q and the diagonal matrix A are the 
skew-symmetry of the former and positive semidefiniteness of the later. No other assumptions have been 
made about a specific form of the matrices Q and A to guarantee the stability of the finite-difference 
scheme [eq. (9)]. 

Remark 3. The discrete operators that are defined by equations [(6)— (8)] are similar in form to those 
that are used for conventional SBP operators (See refs. [9, 10]). What is new, however, is the fact that 
the matrices Q and A depend on u. 

Next, a new one-parameter family of fifth-order WENO schemes is developed, and then is used as the 
starting point in the development of a family of “energy-stable” WENO schemes. Great care is exercised 
in the developing the WENO schemes to ensure that design-order accuracy is achieved in the vicinity of 
smooth extrema. 


3 Fifth- and Sixth-order WENO Schemes 


Any conventional high-order WENO finite-difference scheme for the scalar one-dimensional wave equa- 
tion (1) can be written in the following semidiscrete form: 


du j , fi+i fo-h _ q 
dt Ax 


(13) 


For the fifth-order WENO scheme that is presented in reference [2], the numerical flux / - + 1 is computed 
as a convex combination of three third-order fluxes defined on the following three-point stencils: Sll = 
{xj- 2 , xj-i, Xj}, Sl = {xj~i,Xj,Xj+i}, and Sr = {xj,Xj+i,Xj+ 2 }- (See figure 1.) Note that this set of 
stencils is not symmetric with respect to the ( j + |) point; thus, the fifth-order WENO scheme is biased 


4 




Figure 1: Extended six-point stencil S' 6 , and corresponding candidate stencils Sll, Sl, Sr, and Srr f° r 
one- parameter family of fifth-order WENO schemes. 

in the upwind direction. A central WENO scheme can be constructed from the conventional fifth-order 
WENO scheme by including an additional downwind candidate stencil Srr = {xj+\, Xj+2, Xj+3}, so 
that the collection of all four stencils is symmetric with respect to point (j + |). 1 * * * The WENO flux, 
constructed in this manner, is given by 

fj+h = W J + l/2f^+l/2 + ^+1/2 #+1/2 + <+l/2#bl/2 + <+l/ 2 /#1/2> ( 14 ) 

where /#\/ 2 , r = {LL, L, R, RR} are third-order fluxes defined on these four stencils: 

f{uj- 2 ) \ 


f{uj+ 3 ) / 

and w LL , w L , w R , and w RR are weight functions that are assigned to the four stencils Sll , Sl, Sr, and 
Srr, respectively. The terms w LL , w L , w R , and w RR in equation (14) are nonlinear weight functions. 
These have both preferred values that are derived from an underlying linear scheme as well as solution- 
dependent components. The preferred values are given by 

d LL = jq - ip ; d L = £ - 3p ; d R = ^ + 3p ; d RR = <p , (16) 

where ip is a parameter. The convergence rate of the scheme [eqs. (13)— (16)] with the preferred values 
w M = d ^ and r = {LL, L, R, RR} is greater than or equal to 5 for all values of the parameter ip in 
equation (16). For the specific value p c — < the fifth-order term vanishes and the convergence rate 

1 The above approach to constructing central WENO schemes is proposed in reference [6]. In general, central high-order 

WENO schemes that are built in this manner, are unstable when unresolved features or strong discontinuities exist in the 

computational domain. The generalized energy-stable methodology presented in the next section guarantees the stability 

of the even order approximations while maintaining their nonoscillatory properties. 


-7 11 

-15 2 

2 5 

11 


-1 

-7 2 


f LL ( u j+ 1 / 2 ) \ 
f L ( u j+ 1 / 2 ) 1 

f R {u j+ 1 / 2 ) 6 

f RR (uj+ 1 / 2 ) J 
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increases to 6. The classical fifth-order upwind-biased WENO scheme of Jiang and Shu is obtained for 
cp = 0. 

The weight functions w LL , w L , w R , and w RR needed in equation (14) are given by 


=- a 
o+h E< 


(17) 


where 

a,=d (r) (i + ^|m) , r = {LL,L,R,RR}, (18) 

and e < 0(Aa; 2 ). The functions /W 1 are the classical smoothness indicators. [See equation (59) for the 
general expression.] For fifth-order schemes, they are given by 


( 3 LL = 

if (Wj -2 - 2 u,_i +Ujf 

+ 

\ (Uj - 2 - 4 u,_i + 3 Uj ) 2 

P L = 

13 ( u i— 1 — ^ u j + Uj+ 1) 

+ 

3 (Uj - 1 -Uj+1 f 

P R = 

12 i u j ~~ 2 uj + 1 + Uj+2) 

+ 

| ( 3 Uj - 4 u, + i + Uj +2 ) 2 

pRR = 

if {Uj+i ~ 2u, +2 + u j+ 3 ) 2 

+ 

3 (— 5 u, + i + 8 u j+2 - 3 u, +3 ) 


The t p also varies with discretization order; the expression for T5 is given by 

f (~fj - 2 + 5 /,-i - 10/, + 10/, +1 - 5 f j+2 + /, +3 ) 2 , for ip ± 0 
5 l (fj-2 - 4 /,-i + 6/,- - 4 /, + i + / i+2 ) 2 , for </? = 0 


Expressions for the fourth-, seventh-, and eighth-order WENO schemes appear in the appendix. 

The WENO mechanics expressed in equations (17)-(20) represents a significant departure from the 
mechanics used in the original algorithm by Jiang and Shu [2], For example, equation (18) replaces the 
expression 


d (r) 

ar_ (c + /3W) 2 


( 21 ) 


that is used in the weight functions given in reference [2]. Furthermore, e < 0( Ax 2 ) replaces e = O(10~ 6 ) 
in [2]. Equations (17)-(20) more closely resemble the weight and smoothness indicators proposed by 
Borges et al. in [8] There are differences, however, in how the parameters t p and e are chosen; differences 
motivated by the following observations. 

Borges et al. [8] proposes the following function T5: 


r 5 = \P LL - P R \ 


( 22 ) 


and a fixed value for the parameter e. Although the above weights and smoothness indicators [eqs. 
(17) — (19)] with the value of T5 given by equation (22) significantly outperform the conventional weights 
of Jiang and Shu [2] for both continuous and discontinuous solutions, three remaining shortcomings 
include: 


1. The value of T5 given by equation (22) is not smooth because an absolute value function is used, 
which may reduce accuracy at points where (3 LL — (3 R changes sign. 

2. This approach currently does not generalize to WENO schemes with design orders other than 5. 
For example, neither tq = \/3 LL — /3 R \ nor tq = \f3 LL — j3 RR \ provides the design order of accuracy 
for the central sixth-order WENO scheme. 

3. Near critical points where f, f", and /'" vanish simultaneously, the modified weights [eqs. (17)- 

(19), and (22)] fail to provide the design order of convergence for the fifth-order WENO scheme if 

no constraints are imposed on the parameter e other than e > 0. An example that demonstrates 

this property is presented in Section 6.1. 
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Selecting t$ in equation (20) to be a quadratic function of the fifth-degree undivided difference defined 
on the entire six-point stencil, circumvents the first, and second drawbacks encountered when using the 
scheme suggested in reference [8] . Indeed, T5 is a C°° function in its arguments, and can readily be 
generalized to WENO schemes of order p by choosing t p to be the ptlr-degree undivided difference that 
is defined on the entire (p+ 1)- point stencil. In contrast to the T5 found in reference [8], which is of order 
0(A:r 5 ) for smooth solutions, the proposed function 75 as given by equation (20) is of order 0{ Ax 8 ) and 
0(Aa; 10 ) for ip = 0 and <p =£ 0, respectively; thus much faster convergence of the fifth- and sixth-order 
WENO schemes to the corresponding underlying linear schemes is achieved. 2 

A remedy for the third shortcoming requires an additional modification to the approach proposed 
in reference [8]. Both the fifth- and sixth-order WENO schemes with the new weights that are given 
by equations (17 -20) are design-order accurate for smooth solutions, including points at which the first- 
and second-order derivatives of the solution vanish simultaneously. However, if all derivatives up to 
fourth are equal to zero, then the fifth- and sixth-order WENO schemes locally become only third-order 
accurate. To fully resolve this issue, the constraint e < 0(Ax 2 ) is imposed herein. 

Equations (13)-(20) describe a new family of fifth-order WENO schemes. The primary difference 
between existing WENO schemes and this new family, is in the stencil biasing mechanics described by 
equations (17)-(20). Detailed theoretical justifications for the choice of t p and e used in equations (17)- 
(20) are presented in sections 5.2, 5.3. There, and again in the results section, it is shown that these 
parameters provide fast convergence of the new WENO and ESWENO schemes to the corresponding 
underlying linear schemes for smooth solutions, and deliver improved shock-capturing capabilities near 
unresolved features. 


4 A General Approach to Constructing High-Order Energy- 
Stable Schemes 


Algorithm (1) transforms an existing WENO scheme into an ESWENO scheme that is characterized 
by the following four properties: 1) a bounded energy estimate for arbitrary nonsmooth initial data, 
2) conservation, 3) design order accuracy for sufficiently smooth data, and 4) discontinuity (shock) 
capturing capabilities that are similar to those of the base WENO scheme. 

By construction, algorithm (1) is applicable to 1-D periodic schemes of any order and produces a 
modified scheme that automatically satisfies the first and second properties: bounded energy estimate 
and conservation (See reference [1] for details.) Likewise, as shown in the next section, the additional 
terms added in the ESWENO formulation do not degrade the formal accuracy of the original WENO 
discretization (property three). However, the degree to which the two formulations differ for unresolved 
data is not clear. Thus, the properties of the new ESWENO scheme must be tested to ensure it retains 
the desirable properties of the original formulation. 


ESWENO Schemes of Fifth- and Sixth-order 

We now apply algorithm (1) to the one-parameter family of fifth-order WENO schemes [eqs. (17)— (20)] 
for the scalar 1-D wave equation (1) with a > 0. Combining equation (15) with equation (14) and 
substituting the resulting WENO flux / J+ 1 into equation (13) produces a stencil for the jth grid point 


2 The conventional fifth-order WENO scheme that corresponds to <p = 0, does not include the downwind stencil Srr', 
therefore, the fourth-degree undivided difference built on the available five-point stencil, is included in equation (20). This 
approach to handle the narrow stencil encountered with <p = 0, also generalizes to other orders of accuracy. 
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<) 


Jj 


1 ~ wf - 1/2 ~ <- 1/2 - <- 1/2 - <- 1/2 

2 +< fl /2 + <+ 1/2 + <+ 1/2 + <+ 1/2 


= 0 


(31) 


Because Y^ r w- r ^ = 1, Ag is always equal to zero, and is not included in D\ d . The sign of the diagonal 
terms (A,< i = 1,2,3, could be either positive or negative; thus, the conventional fifth- and sixth- 
order WENO schemes may become locally unstable. If we define (A ®) . ., i = 1,2,3, to be smoothly 
positive 


(A i)jj - otv (Af +$i + 


(Af )jj] 


_ _ . 'J' 

then the additional artificial dissipation operator becomes D^ d = P _1 ^®_i[Di*](Af)[£)i*] , and the 
resulting energy-stable scheme is obtained by adding the additional dissipation term to the original 
WENO scheme. That is, 

D 5 = D 5 + D 5 ad . (32) 


By construction, D 5 satisfies all of the conditions of Theorem 1, thereby providing stability of the fifth- 
and sixth-order ESWENO schemes that are defined by equation (32). 


5 Consistency Analysis 


5.1 Necessary and Sufficient Conditions for Consistency of WENO Schemes 

We now derive the necessary and sufficient conditions for the weight functions u> : for a family of pth- 

order WENO schemes to attain the design order of accuracy. Similar conditions have been obtained 
for the conventional fifth-order WENO scheme in reference [7]. With the same approach discussed in 
section 3 for p = 5, a one-parameter family of pth-order WENO fluxes can be constructed by using a 
convex combination of ( s + 1) fluxes of sth order as 

fj± i/2 = yK+i /2 <h /2 ) (2^) 

r 

with 

/j±i /2 = h(xj± 1 / 2 ) + c i" )Axl + 0(Ax p+1 ), (34) 

l = S 

where h(x) is the numerical flux function that is implicitly defined as 




/ (x) = 


Ax 


Hv) dp 


(35) 


( V ) 

and c\ are constants that do not depend on Ax. 

The corresponding pth-order WENO operator that approximates the first-order spatial derivative is 
given by 


[D p f], = 


f 1+1/2 - fj- 1/2 


E 


(r) 

'j+l/2Jj+l/2 


(r) 


- w (r) f yn 

W j-l/2Jj-l/2) 


(r) 


(36) 


Ax Ax 

where [-]j is a jth component of a vector, the index r in equation (36) sweeps over all (s + 1) stencils, 
and ■uA’A is a nonlinear weight function that is assigned to the corresponding s-point stencil S r . For 
sufficiently smooth solutions, the weights u< approach their preferred values S r \ so that the WENO 
operator converges to the target linear operator U Tar s et as 


-•Target _ .Target E (d^ fj+1/2 ~ d(r) ff} i/ 2 ) 
j-^Targetfj = O+i/2 0_l/ 2 _ r V ' J 


Ax 


Ax 


(37) 
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where dS r > is a one-parameter family of constants chosen to ensure pth-order convergence of the target 
operator to the exact value of the first-order derivative at Xj. That is, 


[_D Target f 


df 


dx 


0{ Ax p ). 


(38) 


The coefficients d ^ for one-parameter families of linear target schemes up to eighth order are given in 
the appendix. All target linear schemes in the family are (2s — l)tli-order accurate, except one central 
scheme, which is (2s)th-order accurate. 

Subtracting equation (37) from equation (36) and using equation (34), we have 


[DP f] . — ^Target fj _ _> 

= S7 E \( w t+i /2 - d(r) )h {r) {x j+1/2 ) - (w[ r } 1/2 - d (r) )/i (r) (* J -_i/ 2 ) 


+ E E Ax ‘ 


J+l/2 
l- lJ r ) 


K'Vl/2 


J~ 1/2 

- dW) - {w { ?}y 2 - d w )l + O(AxP) 


(39) 


From equation (39) it immediately follows that to retain pth-order accuracy, the weights of the WENO 
operator D p should satisfy the following necessary and sufficient conditions: 


E 


W {r) - d M 

W j± 1/2 “ 


= 0{ Ax p+1 ) 


E ^ 1/2 - ^ ) - ( W< f-l/2 ~ ^ ) = °( AxP S+1 ) 

r L J 

E C P ] ( w j r +i/2 - d(r) ) - (^-1/2 - d{r) ) = 0{Ax) 


(40) 


Here, we use the following properties of h(x) and d ^ : f{x) = 1 / 2 ^ ; anc [ ^ d( r ) = 

r 

By construction, the weights [eq. (17)] are normalized such that w- r ' 1 = 1; thus, the first constraint 
in equation (40) is satisfied identically. To simplify the analysis, especially for f(x) with an arbitrary 

(r) 

number of vanishing derivatives, we hereafter use the following sufficient condition on w\ for the 
conventional WENO scheme to attain pth-order accuracy: 

_ d (r) = o(Ax p ~ s+1 ) . (41) 


This constraint is a direct consequence of the necessary and sufficient conditions [eq. (40)] . 


5.2 Sufficient Conditions for Consistency of ESWENO Schemes 

In this section, we show that the conditions [eq. (41)] and the constraints on Si in equation (23): 

Si = 0( Ax p ~ 2i+1 ) , 1 < i < s , (42) 

guarantee that the energy-stable modifications of the conventional pth-order WENO scheme (see section 
4) preserve the design order of the original scheme. As follows from equation (25), the ESWENO operator 
consists of two terms: one is the original WENO operator and the other is the additional artificial 
dissipation operator D a d that is given by equation (24). As shown in the previous section, the WENO 
operator is pth-order accurate if equation (41) holds. Hence, we need only show that D a di = 0( Ax p ). 

Let us prove this conjecture for the one-parameter family of fifth-order ESWENO schemes that is 
presented in Section 4.1. Note that the term P~ 1 D\ A 0 [Hf] T in equation (27) is identically equal to zero 
because of the normalization w\ r ^ = 1; it need not be included in D a d- Thus, the additional artificial 
dissipation term is given by 

D ad f = P~ 1 D 1 A 1 [D 1 ] t { + P- X D\ A 2 [Di] T f + p- 1 D?A 3 [i??] T f (43) 
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i = 1,2,3 



\J \i 2 Sf + A i 


(44) 


where A i is a jth diagonal element of the matrix A.;. For simplicity, we have omitted the superscript 5 
and the subscript (j, j) in this section. 

First, we evaluate the terms Ai, A2, and A3 that are defined by equations (28-30). To simplify the 
derivation, the sufficient conditions [eq. (41)] for the one-parameter family of fifth-order ESWENO 
schemes is rewritten in the following form: 


w[ r) - d (r) 



0(Ax 3 ) + 0( Ax 4 ) . 


(45) 


In equation (45), the stencil width s of each reconstruction polynomial is equal to 3, and the order p of 
this family of schemes is equal to 5 if ip ^ 55 or 6 if As follows from equation (45), the stiffer 

constraint on the weights should be imposed to obtain sixth-order accuracy. 

Replacing w, r ' 1 in A2 with the corresponding preferred values d (r ' ) and using equations (16) and (29), 
for any value of the parameter p yields 

[d LL - 4 d LL + d L -d R + 4 d RR - d RR ] = 0 . (46) 

Subtracting equation (46) from equation (29) and taking into account equation (45) yields 


A2 



LL 

1 + 3/2 


d LL ) - 4(i 


LL 

J j+ 5/2 


d LL ) + (wf +3/ 2 - d ) 


-(<+1/2 - d R ) + 4(w rr 1/2 - d RR ) - ( w rr 1/2 - d RR ) 
(v 5 - m) °(&x 3 ) + 0(Ax 4 ) . 


(47) 


By comparing equations (29) and (30), one can see that Ai is at least one order higher than that of 
A2 because of an additional cancellation that occurs within each group of terms that is associated with 
the same stencil. For example, expanding all of the terms w LL in equation (30) about Xj yields 


^ W j+ 1/2 ^j+3/2 + 2Wj +5 / 2 — 


dw LL 

2 ~lhr 


Ax + 0{ Ax 2 ) 

Xj 


(48) 


which gives an extra factor 0( Ax) compared with the corresponding terms w RL in equation (30): 
w j+ 3/2 ~ dw RR 5 / 2 = 0(1)- The same conclusion can be drawn for the other groups of terms that 
are associated with the L , R 1 and RR stencils in equation (30), which leads to 


Ai 



0(Ax 4 ) + 0(Ax 5 ) . 


(49) 


Applying the same procedure to A3 yields 


A 3 = \ ( d LL ~ d RR ) + \ [(^5/2 - d LL ) - (w rr 1/2 - d RR ) 
= Mlo ~ 2( P) + m) °( Ax3 ) + 0{Ax 4 ) . 


(50) 


To evaluate each term in equation (43), we consider three cases: 1) | A* | » Si > 0, 2) A,; = O(Si), and 3) 
| A* | Si, i— 1,2, 3. If we assume that |A* | » <5, > 0, then equation (44) can be expanded as follows: 


which yields 


\ _ I A: I + A,; | S' 2 

1 _ 2 + 4|Aj| ’ 


Ai = 


Ai 
Si 
<5, 2 

4|Ai| 


(51) 
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if Ai > 0 
if Ai = 0 
if Ai < 0 


(52) 



where the higher order terms have been omitted. Because equation (52) has been derived under the 

8 2 

assumption that |Aj| Si > 0, we can immediately conclude that | A* | Si jjy-j- Therefore, we only 
need consider that 

A* = Xi , (53) 

which provides the lowest order of convergence for the term D\Ki[D\\ T i. If we substitute equation (53) 
in equation (43) and use equations (47, 49, 50), then the additional ESWENO dissipation term becomes 


D a df 

■ 

+ 


Uip-±)0(Ax 3 ) + 0(Ax*) 
[{<p-±)0(Ax 2 ) + 0(Ax 3 ) 


DllDlY f 


To~ 2 ^ 

6Ax 


+ (v - m) 0{Ax 2 ) + 0(Ax 3 ) D 3 [D 3 } T f 


(54) 


If we take into account that D\[D\] T f = 0( Ai 2i ), then D a di can be recast as 


Dadi 



0(Ax 5 ) + 0(Ax 6 ) . 


(55) 


From equation (55) we see, that the additional dissipation term is at least fifth-order accurate for all 
values of the parameter ip. The fifth-order term vanishes for = A- ; for this value the order increase by 
one to sixth-order. 

The second case can be considered in a similar manner. Substituting A, = 0(Si), i = 1,2,3, in 
equation (43) yields 


Dadi 


2m Dl[Dl] T { + Q^ Dl[DV T i+ q^D 3 [D 3 ri 

0(6 iAs) + 0(6 2 Ax 3 ) + 0(S 3 Ax 5 ) . 


(56) 


To guarantee that the ESWENO dissipation term is pth-order accurate, the following constraints should 
be imposed on <5,: 

Sj. = OiAxP- 1 ) 

62 = 0(Ax p ~ 3 ) (57) 

63 = 0(Ax p - 5 ) , 

where p is equal to 5 for ip ^ ^ or 6 for 1 p = Note that the constraints [eq. (57)] are fully consistent 
with those that are given by equation (42). The third case, | A* | <C Si, is similar to the second one and 
results in the same constraints [eq. (57)] on Sp, therefore the third case is not presented here. 

Remark 4. Note that Si, i = 1,2,3, are user-defined parameters; therefore, the conditions [eq. (42)] can 
always be met. 

Remark 5. Although only fifth- and sixth-order ESWENO schemes have been analyzed in this section, 
the same procedure is directly applicable to the additional ESWENO schemes presented in the appendix. 
Thus, we can conclude that if equations (41) and (42) hold, then all the ESWENO schemes that are 
considered in this paper are design-order accurate. 


5.3 Consistency of the ESWENO scheme with New Weights 


The new weight functions for the one-parameter family of pth-order WENO and ESWENO schemes are 
given by 


(V ) 

W K > = 




= d (r) ( 1 + 


-P(r 


(58) 


where /?w) are the classical smoothness indicators: 


S — 1 


/3 (r) =^Aa- 2 '- 1 


1 = 1 


pX ■ 1 1 

f J+% 


' X . 1 

3- -7 


d l q r (x) 

d l x 


dx 


(59) 
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q r (x ) is an (s — l)th-degree reconstruction polynomial that is defined on a stencil S r , and e is a small 
positive parameter that can depend on Ax. In equation (58), t p is defined by 


t p = (V < Xj- s+ i,. . .,x j+s >f , for p 7 0 


(60) 


t p = (V < Xj- S+ 1, . . .,x j+a - 1 >) 2 , for p = 0 


(61) 


where V < Xj- S+ 1, . . . , Xj +S > is the pth-degree, undivided difference. Note that for the original WENO 
schemes of Jiang and Shu [2], which correspond to p = 0, the entire stencil includes only (2s — 1) points; 
therefore, the highest degree of undivided differences that can be constructed on this stencil is 2s — 2, 
rather than 2s — 1, which is used for the other schemes in this one-parameter family. In particular, 
w^ r \ f3^ r \ and 75 for the fifth- and sixth-order WENO and ESWENO schemes are given by equations 
(17)-(20). Another scheme that requires special consideration is the central ESWENO scheme, which is 
obtained by setting p = p c , so that it is one order higher than that of the other schemes in the family. 

First, we present a truncation error analysis for the entire one-parameter family of pth-order ESWENO 
schemes, except for the two schemes that correspond to ip = 0 and p = p c . We demonstrate that the new 
weights that are defined by equations (58)-(60) satisfy the sufficient condition [eq. (41)] for smooth so- 
lutions with any number of vanishing derivatives if the following constraint is imposed on the parameter 
e in equation (58): 


> 0 ( A: 


„p+s— 1 


) > 0, 


(62) 


where p is a design order of the scheme and (s — 1) is a degree of the corresponding reconstruction 
polynomials. 

If we assume that all of the required derivatives are continuous and use the properties of the Newton 
divided differences, then the Taylor series expansions of (3 ^ and t p (tp 7 0) at Xj are given by 


/? (r) = f' 2 Ax 2 + 0(Ax 3 ) , 
Tp = (f ip) ) 2 Ax 2p + 0(Ax 2p+1 ) 


(63) 

(64) 


where f and are the first- and pth-order derivatives of / at Xj. For example, for the family of 
fifth-order WENO schemes, these expansions are 


(3 ll 

(3 L 

p R 

pRR 


= f' 2 Ax 2 + 
= f 2 Ax 2 + 
= f' 2 Ax 2 + 
= f 2 Ax 2 + 


if/" 2 - |/7'") A* 4 + (-f f"f" + \fT) Ax b + 0( Ax 6 ) 


T §/" 2 + I / 7 '") Ax 4 + 0(Ax 6 ) 
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f "2 _ 2 /T „) A;c 4 + (M /T _ Ax 5 + 0(Ax 6 ) 

f" 2 _ Ax 4 + _ 5/'/"") Ax 5 + 0( Ax 6 ) 


(65) 


r 5 = (/ (5) 


Ax 


10 


0( A. 


X 11 ) 


for p 0 . 


( 66 ) 


We first consider a case with f(xj) 7 0. Substituting equations (63) and (64) in equation (58) and 
accounting for equation (62) yields 


./JO) 


= 0( Ax 2p ~ 2 ) (1 - 0(Ax p+s-3 )) , 


(67) 


which leads to 

«,(»•)= d< r ) + 0( A ® 2p “i) . (68) 

Note that the order of convergence of w ^ to its preferred value is 2p — 1 rather than 2 p — 2. The 
main reason for such “superconvergence” is the additional cancellation that occurs because the leading 
truncation error terms of all of the smoothness indicators are identical to each other if f 7 0, as can 
be seen in equation (63). Equations (67) and (68) are valid only for p > 3 and are not applicable to the 
third-order WENO and ESWENO schemes that correspond to p = 0. The detailed analysis of these 
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third-order schemes (p = 0) is presented in reference [1], By comparing equation (68) with equation (41), 
we can immediately conclude that the new weights satisfy the sufficient condition [eq. (41)], ensuring 
that both the WENO and ESWENO schemes are design-order accurate. Furthermore, for schemes of 
order 3 (p ^ 0) or higher, the weights converge to their preferred values at a rate that is significantly 
faster than that given by the sufficient condition (41). The result is a much faster rate of convergence 
for both the ESWENO and WENO schemes to the corresponding target linear schemes, even on coarse 
and moderate grids. 

The next issue that is addressed is the convergence of the ESWENO schemes near the critical points 
at which the first-order and higher order derivatives of the flux approach zero. Let x c be a critical point 
at which the flux function is sufficiently smooth and at which its derivatives of order up to ?i 1)£ jth are 
equal to zero. That is, 

f'(x c) = • • • = f (nvd \x c ) = 0 , f^“ +1 \x c ) ± 0 . 

In contrast to the previous case for which /' ^ 0, the leading truncation error terms of the smoothness 
indicators at the critical point are not equal to each other; thus, no additional cancellation occurs. For 
any number of vanishing derivatives, the following inequalities always hold: 

e > 0{Ax p+s ~ l ) > 0(Ax 2p ) > t p , 


which yields 

t p 

e + /3M 

By using equation (69), the weights can be recast as 

, , d< r > + d (r 

w (r) = 


< 1 . 


i + E 


e +^ (r) = d {r) + o ( Te ^ 


(69) 


(70) 


where we use y^ ir . = 1. For any number of vanishing derivatives, we have 


p ^ 0(Ax 2p J-^ _ 0(Ax p ~ s+1 ) . 


Tp < Tt < 

e + /3(r) - e - 0(Axp+ s ~ 1 ) 


(71) 


From equation (71) we see that w ^ converges to d ^ at the rate of 0(Ax p ~ s+1 ) or higher and satisfies 
the sufficient condition (41). 

As mentioned at the beginning of this section, the two schemes that correspond to p = 0 and p = p c 
require special consideration. Using the same procedure that is outlined above, we can easily show that 
if the parameter e satisfies the constraints 


> 0(Aa; 3s_4 ) , for p = 0 

> 0( Ax 3s ~ 3 ) , for p = p c , 


(72) 


then the corresponding ESWENO schemes are design-order accurate regardless of the number of van- 
ishing derivatives of the solution. The constraints [eq. (72)] are derived using the following relations 
between the order of the scheme and the degree of the reconstruction polynomials: 


p = 2s — 1, for p = 0 
p = 2s, for p = p c . 


Also, note that if no constraint is imposed on e except that it must be strictly positive, then the order 
of convergence of the pth-order WENO and ESWENO schemes with the new weights may deteriorate 
from p to s. Indeed, for a sufficiently large number of vanishing derivatives n v d, (3^ and t p may become 
of the same order. For example, for the family of fifth-order schemes (p ^ 0), this type of degeneration 
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occurs at n v d = 4. If f = f" = f" = f"" = 0, ^ 0, and e < 0(A:r 10 ), which does not satisfy the 

condition given in equation (62), then equations (65) and (66) lead to 


T5 

e + (3 M 


0(A:r 10 ) 

0(Ax 10 ) + 0(Ax 10 ) 


0 ( 1 ) . 


(73) 


As a result, the weights are of order 0(1) and the fifth-order WENO and ESWENO schemes locally 
degenerate to third order. 

Remark 6. In reference [8], the following modification of the weight functions given in equations (17)- 
(19), and (22) is proposed to recover the fifth-order rate of convergence if n v d < 2: 


,,(r) _ 


i 


= d (r 


i + 


T5 


- (30) 


(74) 


where rn = 2. However, the modified weights [eq. (74)] still experience the same degeneration in 
accuracy near critical points for any choice of the parameter m if n v d < 6. We demonstrate this with the 
following “thought” experiment. First, build a polynomial of degree 7 (or higher) such that at least 6 
derivatives vanish at a single point. Next, connect this polynomial to a constant polynomial at the point 
where the derivatives of the first polymonial vanish. The solution to equation (1) is then a six times 
continuously differentiable function everywhere on the combined piecewise polynomials. Note that the 
underlying linear scheme is fifth-order accurate in this case. If we assume that the LL stencil is located 
completely in the constant part of the solution (i.e. fj -2 = fj-i = fj) while the other two stencils 
include points from the polynomial part of the solution, then we have (3 LL = 0, (3 L ^ 0, and (3 R ^ 0. 
As a result, 75 = \(3 LL — (3 R | = (3 R . If no constraint is imposed on e, then we can always choose e such 
that (3 R > £ -> 0 on any given grid, which yields 

This leads to w R = 0(1) and, consequently, to a loss of the design order of accuracy. Note that this 
problem persists regardless of the choice of m in equation (74). 

Remark 7. The parameter e is user-defined, and therefore, the sufficient conditions [eqs. (62) and (72)] 
can always be satisfied. 

Remark 8. Equations (62) and (72) do not provide sharp estimates for the parameter e, which can be 
weakened if additional information regarding the number of vanishing derivatives is available a priori. 
Note, however, that the constraints [eqs. (62) and (72)] guarantee the design order of convergence of 
the corresponding WENO and ESWENO schemes for smooth solutions with an arbitrary number of 
vanishing derivatives. 

The last issue that we discuss in this section concerns discontinuous and unresolved solutions. To 
successfully emulate the ENO strategy, a stencil for which the solution is discontinuous is eliminated 
from the approximation by effectively nullifying the corresponding weight that is associated with this 
stencil. Let a discontinuity be located inside a stencil S r , while the solution is smooth in all other 
stencils. We can easily verify that t p = 0(1), because t p involves all points of the entire stencil including 
those that contain the discontinuity. Therefore, the weight uA r ) can be evaluated as 


w 


(0 




■ 


o(i) + E 

l^r 


0 ( 1 ) 

7+om 

OO) 

e+0(Aj; 2 ) 


0 ( 1 ) 

0(1) + £ +o(aU) 


0(1) [e + 0(Ax 2 ]) 


(75) 


Based on the above eciuation, to reduce the influence of e on the solution near the discontinuity, the 
parameter e should satisfy the following constraint: 

e < 0(Ax 2 ) . (76) 
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Indeed, if equation (76) is met, then w h) is of the order 0( Ax 2 ) and has the same order of magnitude 
as it would have if e = 0. Hence, the parameter e must be bounded not only from below by equations 
(62) and (72), but also from above by equation (76). 

Another consideration that can help us optimally select the parameter e is the manner in which the 
WENO and ESWENO schemes handle small-amplitude oscillations. Consider a solution that contains 
small-amplitude spurious oscillations 

f = fs+Sf , (77) 


where f s is a smooth component and 5f is a nonsmooth, high-frequency component of the solution. By 
substituting equation (77) in equation (59) and assuming that the contribution of the smooth component 
is negligibly small, we have 


/3 (r) = 0(5 f 2 ) . 


As follows from equation (58), if e 7§> 0(5 f 2 ), then the denominator e + 0( r > is dominated by e, and 
these spurious oscillations cannot be detected by the ESWENO dissipation mechanism. However, if 
e < 0(5 f 2 ), then the weights begin to deviate from their preferred values, which increases dissipation in 
those stencils containing spurious oscillations. 

From these considerations it follows that the lower bound of the constraints on e given by equations 
(62), (72), and (76) should be used (1) to obtain the design order of convergence at the smooth parts 
of the solutions, (2) to effectively damp high-frequency spurious oscillations, and (3) to provide good 
shock-capturing capabilities near strong discontinuities. Therefore, in our numerical experiments, we 
use 



0(Ax 3s " 4 ) , 

for p = 0 


0(Ax 3s ~ 3 ) , 

for p = p c 

(78) 

0(Ax 3s - 2 ) , 

otherwise 



where s — 1 is the degree of the reconstruction polynomials. With the above selection of e, all spurious 
oscillations of amplitude 0(e 4 / 2 ) and higher are suppressed by the ESWENO dissipation mechanism. 
The same conclusions can be drawn for the WENO schemes with the new weights given by equations 
(58)-(61). 


6 Numerical Results 


We now assess the performance of the fourth-, fifth-, and sixth-order ESWENO schemes with the new 
weights and compare them with the conventional WENO counterparts. For all of the numerical experi- 
ments that are presented, the parameters e of the ESWENO weight functions and the parameters <5, of 
the additional artificial dissipation operator are chosen based on equations (78) and (42) with two minor 
modifications. First, to preserve the scale invariance of the original WENO schemes, the physical grid 
spacing Ax in equations (42) and (78) is replaced with A£, where A£ = 1/J and J is the total number of 
grid cells. Second, as follows from equation (58), the parameter e should be scaled consistently with pO\ 
In regions where the solution is smooth, for the third- and fourth-order WENO and ESWENO schemes, 
/ 3 ( r ' ) approximates j' 2 A£ 2 . For the fifth- and sixth-order schemes, (3 ^ approximates the linear combi- 
nation of the same first-order derivative term and a second-order derivative term that is proportional to 
f" 2 . The same pattern persists for higher order schemes as well. In the vicinity of discontinuities, P^ 
is of the order 0(f 2 ). If we take into account the above considerations, for all test problems considered 
the parameter e is set to 


e = 
C 


r ca£ 3s - 4 , 

for ip = 0 


\ CA£ 3s ~ 3 , 

for p = p c 


[ CA£ 3s-2 , 

otherwise 

(79) 

max (||/o II- ll/o 2 ||, - - ■ , II (fo 1} ) II j . 



where /o = /[«o(£)] is the initial flux, uq(£) is the initial condition, /q S ^ is the (s— l)th-order derivative 
of /o with respect to £, s — 1 is the degree of the reconstruction polynomials, pj is a set of points at 
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exact 



Figure 2: Solutions obtained with fourth-order WENO and ESWENO schemes on 101-point grid for 
linear wave equation with initial condition [eq. (81)] at t = 0.04. 

which the solution is discontinuous, and || • || is a norm in which the solution is sought. The scaling factor 
C in equation (79) can easily be evaluated because it depends only on the initial condition for which 
the locations of all discontinuities are known a priori. Note that the parameter e is calculated once and 
the same value is used over the entire time interval of integration; thus the computational cost does not 
increase. 

In accordance with equation (23), the parameter <5, should be scaled consistently with A *. Because 
A i is a linear combination of the weights with each being of the order 0(1), the scaling factor of Si is set 
equal to one, which leads to 

Si = A fP- 2i+1 . (80) 

Equations (79) and (80) eliminate the ambiguity in determining the parameters e and S, for the ESWENO 
schemes; thus, the equations for e and S t are free of tuning parameters. Furthermore, equations (79) 
and (80) are fully consistent with the sufficient conditions [eqs. (42) and (78)] and allow the ESWENO 
schemes to be invariant when the spatial and time variables are scaled by the same factor. Note that 
the parameter e for the conventional WENO schemes of Jiang and Shu is set to 10~ 6 as recommended 
in [2]. In accordance with the recommendations of Borges et al. [8], the parameter e for the fifth-order 
WENO scheme with the weights given by equations (17) — (19), and (22), which is referred as WENO-Z, 
is set to 10“ 40 . 

The time derivative for all steady test problems is approximated by using a third-order total variation 
diminishing (TVD) Runge-Kutta method that is developed in reference [11], while unsteady problems 
are integrated by using a fourth-order low-storage Runge-Kutta method (ref. [12]). To reduce the fourth- 
order temporal error component and make it consistent with the spatial error of fifth- and sixth-order 
schemes, the time step in global grid refinement studies is reduced by a factor of 2 6 / 4 for each doubling 
of the number of grid points in space. The Courant-Friedrich-Levy (CFL) number has been set to 0.3 
and 0.6 for the steady and unsteady test problems, respectively. 


6.1 Scalar Linear Wave Equation 

We begin by verifying that the new class of ESWENO schemes provide the design order of convergence 
for smooth problems, including local extrema. To check this property, we consider equation (1) with 
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a = 1 and the following initial condition: 


Uo (x) = e- 300 (*-*«) 3 , (81) 

where x c is 0.5. The computational domain for this test problem is set to 0 < x < 1. Numerical solutions 
are calculated on a sequence of globally refined uniform grids and advanced in time up to f = 1, which 
corresponds to one period in time. 

First, we show that the conventional fourth-order WENO scheme is unstable for this smooth problem; 
however, the corresponding ESWENO is stable, and its solution is in excellent agreement with the 
exact solution, as shown in Figure 2. To make the conventional fourth-order WENO scheme stable, a 



Figure 3: L 0 0 error norms obtained with the fourth-order linear, WENO, and ESWENO schemes for 
linear wave equation with initial condition given in equation (81). 


smoothness indicator that corresponds to the downwind stencil has been modified as follows: 


(3 R = 


(, 0 L ) k + (f3 c ) k + ( (3 R ) k 


1 1/fc 


(82) 


where k is a constant that is greater than 1. Note that j3 R is a C°° function of its arguments and 
approaches max(/3 L , (3 C , (3 R ) as k — ► oo. In contrast to equation (82), a modified smoothness indicator 
j3 R = max(/3 L , (3 C , f3 R ) that is proposed in reference [6] is a nonsmooth function of (3 L , (3 C and /3 fl , which 
may lead to the degeneration of the design order of accuracy and is more prone to spurious oscillations 
near unresolved features and strong discontinuities. For k — *■ oo, the downwind smoothness indicator 
that is given by equation (82) prevents the corresponding weight w R from being larger than the lesser 
of the other two weight functions; thus, the stencil is biased in the upwind direction near unresolved 
features. In smooth regions, all three smoothness indicators are of the same order, and the weights 
approach their preferred values of w L = g, w c = |, and w R = g, which provides the design order of 
accuracy. For all of the numerical experiments that are presented herein, the parameter k in equation 
(82) is equal to 4. 

Figure 3 shows the L error norms that are obtained with the fourth-order WENO and ESWENO 
schemes, and the corresponding underlying linear schemes. As shown in Figure 3, the fourth-order 
ESWENO scheme is significantly more accurate than the conventional fourth-order WENO scheme. The 
Loo error that is obtained with the fourth-order ESWENO scheme is slightly higher than that obtained 
with the corresponding linear scheme on coarse meshes and reaches its theoretical limit starting at 
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Figure 4: error norms obtained with the fifth- (left) and sixth-order (right) linear, WENO, WENO-Z, 

and ESWENO schemes for linear wave equation with initial condition given in equation (82). 


J = 201 grid points. The maximum error occurs at the peak of the Gaussian pulse indicating that the 
ESWENO scheme is design-order accurate at the smooth extrema. In contrast to the ESWENO scheme, 
the conventional WENO scheme demonstrates only a third-order convergence rate, even on the finest 
mesh with J = 1601, and is two to three orders of magnitude less accurate on moderate and fine meshes. 

If the new weights [eqs. (58) and (59)] are used instead of their conventional counterparts, then the 
design order of convergence of the fourth-order central WENO scheme is recovered, as shown in figure 
3. Moreover, the fourth-order WENO scheme with the new weights provides slightly better accuracy 
on coarse meshes than the corresponding ESWENO scheme. This result is not surprising because the 
ESWENO scheme has the additional dissipation term [eq. (24)], which guarantees the stability. In 
general, if a WENO scheme with the new weights is stable, then it is slightly less dissipative than the 
corresponding ESWENO counterpart. For all of the problems that are considered, however, the results 
obtained with the conventional WENO schemes with the weights given by equations (58) and (59) are 
practically indistinguishable from those of the corresponding ESWENO schemes; therefore these results 
are not presented hereafter. 

The Lqo error norms that are obtained with the fifth-order WENO-Z scheme and the fifth- and 
sixth-order WENO and ESWENO schemes for the same test problem are depicted in figure 4. Similar 
to the fourth-order case, the conventional sixth-order central WENO scheme is unstable. This problem 
is avoided by using the same upwinding technique that is outlined earlier. For the sixth-order central 
WENO scheme, the modified smoothness indicator that corresponds to the most downwind stencil is 
given by 


/ 3 rr = 


{/3 LL ) k + ( (3 L ) k + (, g R ) k + ( [3 RR ) k 


1 1/fe 


(83) 


Both the fifth-order WENO-Z and fifth- and sixth-order ESWENO schemes are equal in accuracy 
to the corresponding underlying linear schemes, as shown in Figure 4. Although the fifth-order WENO 
scheme exhibits the design-order convergence rate, its L 0 0 error norm is nearly an order of magnitude 
larger than those of the corresponding fifth-order WENO-Z and ESWENO schemes. In contrast to that 
of the sixth-order ESWENO scheme, the L 0 0 error that is obtained with the conventional sixth-order 
WENO scheme shows only fifth-order convergence and is quite far from the theoretical limit (represented 
by the corresponding sixth-order central linear scheme). Note that the WENO-Z schemes of fourth- and 
sixth-order are currently unavailable in the published literature and therefore are not presented herein. 

As we have proven in sections 2 and 4.1, for the linear convection equation (1) with piecewise 
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5th-order WENO 
5th-order ESWENO 



Figure 5: Time histories of rightmost eigenvalue of symmetric part of fifth-order WENO and ESWENO 
operators for Gaussian pulse problem. 


continuous initial conditions, the family of ESWENO schemes is stable in the energy norm. The sufficient 
condition for stability is the negative semidefiniteness of — ^(D + D 1 "), where D is the ESWENO discrete 
operator. This property implies that all eigenvalues of the symmetric portion of the ESWENO operator 
are nonpositive. In contradistinction to the ESWENO scheme, the symmetric portion of the WENO 
operator may have positive eigenvalues. (See section 4.1.) Note that the same conclusion can be drawn 
for the fifth-order WENO-Z scheme (ref. [8]), whose derivative operator is identical to that of the 
conventional WENO scheme with the modified weights [eqs. (17) — (19), and (22)] which vary in the 
interval from 0 to 1 near the unresolved features. These properties of the WENO, WENO-Z, and 
ESWENO schemes are shown in figure 5, which gives the time histories of the rightmost eigenvalue 
of the symmetric portion of the operators computed on a 201-point grid. The rightmost eigenvalue of 
the ESWENO operator — ^(D + D T ) is equal to zero up to the order of the round-off error, while the 
symmetric part of the conventional fifth-order WENO and WENO-Z operators have positive eigenvalues 
of order (^(lO” 1 ) and 0(10), respectively. (See figure 5.) This results from the fact that the symmetric 
portion of these WENO-type operators is not negative semidefinite, if unresolved features exist in the 
computational domain. Note that the presence of the positive eigenvalues does not imply that the 
fifth-order WENO and WENO-Z schemes are globally unstable because these eigenvalues correspond to 
different grid points at different moments of time. 

The superiority of the ESWENO schemes with the new weights compared with the conventional and 
modified WENO counterparts, becomes evident when a more challenging problem is considered. In the 
previous test problem, the initial condition [eq. (81)] has the critical point at which f = 0 and /" ^ 0 
(i.e. , the number of vanishing derivatives is n v d = 1). As shown in section 5.3, the fifth-order WENO-Z 
scheme that was developed by Borges et al. [8] fails to recover the design order of convergence if n v d > 6 
for an arbitrary choice of the parameter m in equation (74). We now demonstrate this property for 
Ti v d > 3. Consider equation (1) with the following initial condition: 

r z 18 - 14z 16 + 69z 14 - 175z 12 + 259z 10 

u 0 (x) = < — 231z 8 + 119z 6 — 29z 4 + 1, for |z| < 1 (84) 

[ 0 otherwise 

where z = 5(x — 0.5) and 0 < x < 1. The above function is six times continuously differentiable 
and has three critical points: one at x = 0.5 with n v d = 3 (i.e., /'( 0.5) = /"(0.5) = /'"(0.5) = 0 and 
f""{ 0.5) ^ 0) and the remaining two at x = 0.3 and 0.7 with n v d = 6. If we compare the L 0 0 error norms 
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Figure 6: L 0 0 error norms obtained with fifth-order linear, WENO, WENO-Z, and ESWENO schemes 
for linear wave equation with initial conditions given in equation (84) at t = 1. 


computed with the fifth-order linear, WENO, WENO-Z, and ESWENO schemes at t = 1, (see fig. 6), we 
see that the situation changes dramatically as compared with the previous test problem. As expected, 
the WENO-Z scheme fails to deliver fifth order convergence, while the ESWENO scheme is fifth-order 
accurate and provides practically the same error convergence as the underlying linear scheme. This 
result is not surprising because the parameter e for the WENO-Z scheme is set to 10 -40 , as suggested 
in reference [8]. As a result, at a point x = 0.5 + Ax the weights become of order 0(1). Thus, if we use 
equations [(65), (22), and (74)] with e = 0 and take into account that for the polynomial in equation 
(84), f = 0( Ax 3 ), f" = 0(Ax 2 ), f" = 0( Ax), and /"" = 0(1) at x = 0.5 + Ax, then we have 


w ( r ) _ d (r) = q 


\ff"f'"~ff""\Ax 5 
f' 2 Ax 2 + (f f' 2 - /'/"') Aa 



(85) 


As follows from the above estimate, the WENO-Z scheme fails to recover fifth order convergence regard- 
less of the choice of the parameter m in equation (74). Our numerical results corroborate this conclusion. 
Figure 6 shows that increasing the exponent in in equation (74) results in an even larger deterioration in 
accuracy. In contrast to the WENO-Z scheme, the conventional fifth-order scheme of Jiang and Shu[2] 
converges at the design order and begins to approach the theoretical limit on the finest meshes. The 
main reason for such a behavior is the choice of the parameter e, which is set to 10 “ 6 . As the grid is 
refined, (3 ^ — > 0, the denominator in equation (21) is dominated by e, and the weights approach their 
preferred values. 


6.2 The 1-D Euler Equations 

In reference [1], the third-order ESWENO scheme is proved to be energy-stable for a system of linear 
hyperbolic equations with periodic boundary conditions. This result can be directly extended to the 
class of higher order ESWENO schemes that are presented in this paper. For nonlinear conservation 
laws with nonperiodic boundary conditions, a similar proof for the high-order ESWENO schemes is not 
currently available. Nevertheless, we would like to test how the new schemes perform for the system of 
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the quasi- 1-D Euler equations, which are given by 





<9U ■ d¥ . 
dt ' dx _ 

= G 
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pu z 

_ (. E + P)u _ 


P= (7-1) (£+*£) , 


( 86 ) 


where A = A(x) is the cross-sectional area of a quasi-l-D nozzle and 7 = 1.4. As has been shown 
in reference [1], the ESWENO reconstruction must be implemented in local characteristic fields to 
guarantee the stability of the ESWENO scheme for systems of hyperbolic conservation laws. In the 
present analysis, both the WENO and ESWENO reconstructions are based on the Lax-Friedrichs flux 
splitting. See reference [1] for further details. 


3rd-order linear 



5th-order linear 



iog,„(N) 


Figure 7: error norms obtained with the third- and fourth-order and fifth- and sixth-order linear, 

WENO, and ESWENO schemes for the subsonic quasi-l-D nozzle problem. 

For the 1-D Euler equations, the preferred biasing of the stencil in the upwind direction, given by 
equations (82) and (83), is used for the central fourth- and sixth-order WENO and ESWENO schemes 
to suppress spurious oscillations near strong discontinuities. In the ESWENO case, these oscillations can 
also be eliminated by increasing the coefficient in front of the D±AiDf term in the artificial dissipation 
operator. However, this approach is more dissipative and also reduces the maximum CFL number for 
which the scheme remains stable. For the test problems that are presented in this section, the results 
obtained with the fifth-order WENO-Z scheme are practically indistinguishable from those of the fifth- 
order ESWENO scheme; therefore the WENO-Z results are not presented hereafter. 

To verify that the new ESWENO schemes are design-order accurate for hyperbolic systems, we 
consider the steady-state isentropic flow through a quasi-l-D nozzle with the following cross-sectional 
area: A(x) = 1 — 0.8x(l — x), 0 < x < 1, as a test problem. The inflow Mach number is set to 0.5 and 
the pressure at x = 1 is assumed to be equal to that at x = 0. Under these conditions, the flow is fully 
subsonic, and the solution is smooth. Global grid-refinement studies for the third-, fourth-, fifth-, and 
sixth-order WENO and ESWENO schemes are presented in figure 7. 

The Loo error norms that are obtained with the third- and fourth-order ESWENO schemes exhibit 
the design order of convergence. On the finest mesh with J = 801, the fourth-order ESWENO solution is 
approximately three and four orders of magnitude more accurate than those of the third-order ESWENO 
and WENO schemes, respectively. As shown in figure 7, the fifth-order ESWENO scheme provides the 
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same accuracy as its underlying linear scheme, thereby exhibiting the perfect error convergence for 
this test problem. Although the conventional fifth-order WENO scheme exhibits the design order of 
convergence, its error norm is larger by a factor of 3 than that obtained with ESWENO counterpart. 
Similar to the fourth-order case, the central sixth-order ESWENO scheme converges at the design-order 
rate on both coarse and fine meshes, thus providing significantly more accurate solutions than both 
fifth-order schemes. 

Remark 9. For this problem, the conventional central fourth- and sixth-order WENO schemes do not 
converge to a steady-state solution even with the upwinding mechanisms [eqs. (82) and (83)] turned 
on. The primary reason for this tendency is the presence of positive eigenvalues in the spectrum of the 
symmetric part of the conventional WENO operator. 



Figure 8: Comparison of fifth-order WENO and ESWENO schemes for steady transonic flow through 
quasi-l-D nozzle. 

The next test problem is the steady transonic flow through a quasi-l-D nozzle with the following 
cross-sectional area: 

A{x) = 1.398 + 0.347 tanh(0.8x — 4), 0 < x < 10. 

The Mach number at x = 0 is 1.5, and the outflow conditions have been chosen so that the shock is 
located at x = 5. Density profiles are calculated using the fifth-order WENO and fifth- and sixth-order 
ESWENO schemes on a 51-point grid and are compared with the exact solution in figure 8. For all 
schemes, the captured shock is smeared over two grid cells, and the numerical solutions are essentially 
nonoscillatory and agree quite well with the exact solution. Note that the fifth- and sixth-order ESWENO 
schemes converge to the steady-state solution an order of magnitude faster than the conventional fifth- 
order WENO scheme; this result again indicates the presence of unstable modes that are generated by 
the positive eigenvalues of the WENO dissipation operator. 

The last two problems considered are standard unsteady problems for testing shock-capturing schemes. 
The first problem is the Riemann problem with the initial conditions proposed by Sod[13]: 

, . f (1,0,1) if -0.5 < x < 0 

(p, u, f) - j ^ 0 _ 125) Q) 0 if 0 < x < 0.5 

The numerical solutions that are computed with both fourth- and sixth-order WENO and ESWENO 
schemes are essentially nonoscillatory, as one can see in figure 9. Note, however, that the fourth- and 
sixth-order ESWENO schemes provide better resolution near the contact discontinuity and the shock as 
compared with the conventional WENO schemes. 
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Figure 9: Density profiles computed with the fourth-order (left) and sixth-order (right) WENO and 
ESWENO schemes on a uniform grid with 201 points for the Sod problem at t = 0.16. 


We conclude this section by presenting the numerical results that are obtained with WENO and 
ESWENO schemes for the shock entropy- wave interaction problem. The solution of this benchmark 
problem contains both strong discontinuities and smooth structures, and is well suited for testing high- 
order shock-capturing schemes. The governing equations are the time-dependent 1-D Euler equations 
[eq. (86)] with G = 0, subject to the following initial conditions: 


( (3.857134,2.629369,10.33333), if - 5 < x < -4 

\ (1 + 0.2 sin 5x, 0,1), if - 4 < x < 5 . 


(87) 


The governing equations are integrated in time up to t = 1.8. The exact solution to this problem is 
not available. Therefore, a numerical solution that is obtained with the conventional fifth-order WENO 
scheme on a uniform grid with J = 4001 grid points is used as a reference solution. 

Numerical solutions that are computed with the fourth-, fifth-, and sixth-order WENO and ESWENO 
schemes on a 301-point grid at t = 1.8 are compared with the “exact” reference solution in figure 10. 
Both the WENO and ESWENO solutions are free of spurious oscillations. The fourth-order WENO 
scheme is the most dissipative among the considered schemes and has the worst resolution in a region 
just upstream of the moving shock. The fourth-order ESWENO scheme provides better resolution of the 
high-frequency oscillations behind the shock. However, the wave amplitudes are much lower than those 
that are computed with the fifth- and sixth-order schemes. Again, the fifth- and sixth-order ESWENO 
solutions are more accurate than those of the corresponding WENO schemes. Both ESWENO schemes 
resolve the smooth extrema quite well, and the solutions are essentially nonoscillatory near the captured 
shock. 


7 Conclusions 

A systematic methodology is developed to construct Energy Stable weighted essentially nonoscillatory 
(ESWENO) finite-difference schemes of arbitrary order, and is used to construct ESWENO schemes 
based on the WENO schemes that are presented in references [2] and [6]. The new ESWENO schemes 
differ from the conventional WENO schemes by the addition of a nonlinear artificial dissipation term 
of special form. The additional term is design-order accurate for smooth solutions, including smooth 
extrema, and guarantees that the ESWENO scheme is stable in the .^-energy norm for both continu- 
ous and discontinuous solutions of hyperbolic systems; for the conventional WENO scheme, an energy 
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exact exact 




6th-order WENO 



Figure 10: Density profiles computed with the fourth-order (top left), fifth-order (top right), and sixth- 
order (bottom) WENO and ESWENO schemes on uniform grid with 301 grid points for density sine 
wave/shock- interaction problem. 


estimate is not available. The distinctive feature of the new class of ESWENO schemes is that the spec- 
trum of the symmetric part of the ESWENO operator is always located in the left half-plane, while the 
symmetric part of the WENO operator has positive eigenvalues; thus, the conventional WENO schemes 
may become locally unstable. 


New weight functions are also developed, guided by newly derived sufficient conditions guaranteeing 
design order accuracy for arbitrary WENO schemes; sufficient conditions that are valid regardless of 
the number of vanishing solution derivatives. A truncation error analysis is used to optimize the weight 
functions, so that the new ESWENO schemes satisfy the sufficient conditions for accuracy and provide 
essentially nonoscillatory solutions for problems with strong discontinuities. Numerical experiments 
confirm that the high-order ESWENO schemes with the new weights are significantly more accurate than 
conventional WENO schemes of the same order, and demonstrate excellent shock-capturing capabilities. 
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8 Appendix 

8.1 Existence of Symmetric Decomposition 

• T 

The following lemma establishes the existence of the decomposition D sym = XAoI-^i^Ad-Di*] • The 
lemma provides an algorithm to decompose a matrix of arbitrary bandwidth 2s + 1 into the sum of two 
matrices: one that can be expressed as one term of the desired symmetric form, and one of bandwidth 
2 (s — 1) + 1. 

Lemma 1. Define E s to be an TV-dimensional, symmetric, bidiagonal matrix with the subdiagonal ele- 
ments e s + yi, e s _)_2,2j • • • , ejv-pjv-s-i, £n,n-s- The parameter s is the offset from the main diagonal. 
Further, define B s _ \ to be an TV-dimensional, symmetric, banded matrix of bandwidth 2(s — 1) + 1 
with elements that are functions of the matrix E s . The matrices E s and B s _ \ satisfy the following 
relationship: 

E s = (-1 ) s [ZV]A s [TV] T + .B s _ 1 . 

with A s = diag[e s +i^i, e s +2,2, • • • , ejv-i,iv-s-i, £n,n-s, 0 s t ] and 0 S is a zero vector of dimension s. 

Proof: A proof by induction is not included herein. One step of the algorithm can be verified 
immediately by inspection. 

A simple example illustrates Lemma 1. Define the 5x5 matrices E 2 , D 15, and B\ such that 


= (-1) 2 [^15 
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Remark 10. Lemma 1 can be used recursively, beginning with the outermost s diagonal and working 
inward, to establish that D sym = 5Z* =0 [-Di*]Aj[.Di®] . Each step generates a new symmetric matrix 
with a bandwidth that is equal to 2 less than the previous value until only a diagonal matrix remains. 
Thus, the algorithm terminates after s + 1 steps with the desired decomposition. 

Remark 11. This recursive algorithm works for nonperiodic matrices in multiple spatial dimensions. 


8.2 ESWENO Schemes of Third- and Fourth-Order 

Based on equation (13), fj + 1 defined for a third- or fourth-order WENO scheme is 



- W j+l/ 2 fj+l /2 + W j+l/ 2 fj+l /2 + W j+l/ 2 fj+l /2 


( 88 ) 


where w L , w c , and w R are weight functions that are assigned to each respective stencil. The second- 
order linear fluxes: /’y+ 1 / 2 , defined in equation (88) for r = { L , C, R}, are 


f L ( u j+ 1/2) 
f C ( u j+ 1/2) 

,f R (u j+ 1/2) 


/ -1 3 0 0 

0 110 
\ 0 0 3 -1 


/ f(uj- 1) \ 
f( u j ) 
f(Uj+ 1) 

V f(u j+ 2) ) 


(89) 
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The terms w L , w c , and w R are the nonlinear weight functions that are given by equations (58) and 
(59) with 

_ f {-fj-i +3/j -3/j+i +/i+ 2 ) 2 , for^yfO 

3 l if 3 - 1 - 2 /j + fj+ 1) 2 . for ip = 0 

(See ref. [1] for further details). These weight functions have preferred values that are derived from 
underlying linear schemes, as well as solution-dependent components. The preferred values are given by 
the formulae 

df = ^ — cp ; dp = | ; d R = <p , (90) 

where <p is a parameter. The convergence rate of equation (13), with the preferred weight values that 
are given in equation (88), is equal to 3 for all values of the parameter ip except p> c = for which the 
convergence rate is 4. 

Combining equation (89) with equation (88) and substituting the resulting WENO flux f J+ i into 
equation (13) produces a stencil for the j th grid point of the form 


D h 
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2Ax 
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V 1/2 
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~“P+ 1/2 
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- 3w f-l/2 
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3w f+l/2 

W f+l/2 
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-<- 1/2 
- 3 <— 1/2 


0 

«2 

3 <+l/2 

0 
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3 - 1/2 


0 

0 

-<+1/2 

0 

0 

0 


( f(uj- 2 ) \ 
f(uj- 1 ) 
f(uj ) 
/K'+i) 

V /(«j+2) / 


(91) 


with fc on the interval j — 2 le k < j+2. 
Simplifying Zl 3 using the expression 


w. 


= 1 - w L - 


yields expressions for a single row of the D 3 matrix of the following form 


D 3 

u j,k 


2Ax 


w 


L 2 <-y 2 


W 

- up 


J — 1/2 




1 


P + 1/2 ■+■ — 1/2 

p — 1/2 2u, j+l/2 ~ 2 W j-l/2 ~ W j+ 1/2 

■< +1/ 2 + <-1^2 + 2 <+l/2 + 1 


U T+l/2 


(92) 


The derivative matrix Z? 3 is now decomposed into symmetric and skew-symmetric parts 

n 3 = n 3 + n 3 

skew 1 sym 

As with the fifth-order case, the skew-symmetric component and the norm of D 3 take the forms 


n 3 

skew 


= P-'Qs ; Qs + Ql = 0 ; P = Axl . 


while the matrix D 3 is expressed as 


D 3 = P- 1 

sym 


( D 1 2 A 3 [D^f + D x x A 3 [.D/f + D x ° A 3 [D^f ) 


(93) 
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where is a diagonal matrix with the expressions for the jth element defined by 



^2 - 4 <+3/2 *4+1/2 

(94) 

0+7 

1 

” 4 

<+3/2 + <+1/2 <+1/2 + W f- 1/2 

(95) 

( X 0 

1 

2 

[ -<_ 1/2 - < 1/2 - <- 1/2 

(96) 



+<+ 1/2 + <+ 1/2 + <+ 1 / 2 j = 0 • 


If we define (Af ) . i = 1,2, to be smoothly positive, then 

= ^[V / (A?)L+^ 2 - , 

and the additional artificial dissipation operator is given by 

2 

i=l 

The resulting energy-stable scheme is obtained by adding the additional dissipation term to the corre- 
sponding conventional WENO scheme as 


D 3 = 


D 3 


D 


ad 


8.3 ESWENO Schemes of Seventh- and Eighth-order 

The fj + 1 , defined for the seventh- and eighth-order WENO schemes, is 


f j+ i 


+ 


w 

w 


LL fLL 
j+l/ 2 Jj+l /2 
R fR 
j+l/ 2 Jj+l /2 


(97) 


where w LL , w L , w c , w R , and w RR are weight functions that are assigned to each respective stencil, 
and are given by equations (58) and (59) with 


f (/,_ 3 - 7 /,_ 2 + 21/, _i - 35 fj + 35 f j+1 - 21 f j+2 + 7 f j+3 - /,+ 4 ) 2 , for <p ± 0 
\ (~fj - 3 + Gfj - 2 - 15 /, -1 + 20/, - 15/, +1 + 6/, +2 - /, +3 ) 2 , for y = 0 


The fourth-order linear fluxes « f° r r = {LL, L,C, R, RR}, which are used in equation (97) are 
defined as 
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( 98 ) 


The derivative matrix D' is now decomposed into symmetric and skew-symmetric parts as 


D 7 


= D 


skew 


+ D 


7 

sym • 
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As with the fifth-order case, the skew-symmetric component of D 7 takes the form 


D 7 

^ skew 


= P~ 1 Qt, Qt + Qj= 0; P = Axl . 

The scheme converges with seventh-order accuracy if the weights are equal to their preferred values 


d^ L = = §— 8 <p,d C = i ,d R = ± + 8<p,d RR = 

and with eighth-order accuracy for the specific value ip c = 

The D 7 sym can be simplified into the form 


( 99 ) 
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where A 7 is a diagonal matrix with expressions for the jth element defined by 


= 


,„ LL _ ,,,RR 

W j+ 7/2 W j+ 1/2 


(^)jj 24 


+ w j+5/2 9 W j+ 7/2 + W j+ 5/2 


_W j+l/2 + /■> - 


4-1/2 


■u; 


rr 


4+1/2 


= Si + 2w f+3/2 “ llu, f-W5/2 + 9u> : 


4+3/2 
+2'ui i 


^ 4 + 5/2 

-2w L 


LL 

°j+ 7/2 


~ 24 


4+3/2 “ ^4+5/2 
- W 4 C +l/2 + U f+3/2 

+ 2u f- 1/2 - 2w f+l/2 
-9<4/2 + H<-l/2 - 2 <+l/2 


+ 6w f+ Z i/2 - 13u #3/2 + 10w f+ i 5/2 - 3 w #7/2 

4” 3 w; *+1/2 — ^1+3/2 + W i+5/2 
' " U ?-l/2 ~ W i+3/2 


_W 4-3/2 + 4 ^-i /o - 3 «4 


i- 1/2 


A+l/2 


+3^4/2 - 10u 41/2 + 13 ^-V2 - 6w^ /2 


(*0)4,4 = 3 [ - U; f-l /2 - <4-1/2 - <4-1/2 - <4-1/2 - <44/2 


+ W i+ 1/2 + <4-i /o + w£_i /o + 41 1 /o + W. 


n L 

L i+ 1/2 


A+ 1/2 


A+ 1/2 


RR 

*+1/2 


= 0 


If we Define (A?) • ,, i = 1, 4, to be smoothly positive as 


( 100 ) 


( 101 ) 

( 102 ) 

( 103 ) 


( 104 ) 


( 105 ) 


(Aj )j,j = -WWrjj+St - W)fj], 


the additional artificial dissipation operator is given by 

4 


D 7 ad = P-^D^^iDS] 


i=l 


and the resulting energy-stable scheme is obtained by adding the additional dissipation term to the 
corresponding conventional WENO scheme as 


D 7 = D 7 


D 7 ad ■ 
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